Uniform

Here is an example of using projectq2a to sample from the standard uniform distribution.

The uniform distribution has probability density function:

\[\text{pdf = }\begin{cases}\frac{1}{b - a} \text{for } x \in [a, b]\\0 \text{ otherwise}\end{cases}\] We first generate a bunch of deviates and then plot their density in teal Then the true density is overlayed as a yellow line.

df_uniform_samples <- 
  data.frame(u = projectq2a(n = 10000, 
                       pdf = dunif, 
                       a = 0, 
                       b = 1, 
                       C = 1, 
                       min = 0, 
                       max = 1))

ggplot(df_uniform_samples, aes(x = u)) + 
  geom_density(color = NA, fill = "#35a09c") + 
  stat_function(aes(x = 0), fun = dunif, args = list(min = 0, max = 1), color = "#fff735", size = 1) + 
  labs(x = "X", 
       y = "Density") + 
  theme_minimal()

Beta

Here is an example of using projectq2a to sample from a beta distribution with both shape parameters equal to 2.

The beta distribution has probability density function:

\[\text{pdf = }\begin{cases}\frac{x^{\alpha - 1}(1 - x)^{\beta - 1}}{B(\alpha, \beta)} \text{ for }0 \le x \le 1\\0 \text{ otherwise}\end{cases}\]

\[\text{where }B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}\text{ with }\alpha, \beta \gt 0\]

We first generate a bunch of deviates and then plot their density in teal. Then the true density is overlayed as a yellow line.

df_beta_samples <- 
  data.frame(b = projectq2a(n = 10000, 
                       pdf = dbeta, 
                       a = 0, 
                       b = 1, 
                       C = 1.5, 
                       shape1 = 2, 
                       shape2 = 2))

ggplot(df_beta_samples, aes(x = b)) + 
  geom_density(color = NA, fill = "#35a09c") + 
  stat_function(aes(x = 0), fun = dbeta, args = list(shape1 = 2, shape2 = 2), color = "#fff735", size = 1) + 
  labs(x = "X", 
       y = "Density") + 
  theme_minimal()

Custom

Here is an example of using projectq2a to sample from a custom distribution.

In this case, the custom density function is that of a beta with shape parameters 2 and 5.

\[\text{pdf = }\begin{cases}30x(1 - x)^4 \text{ for }0 \le x \le 1\\0 \text{ otherwise}\end{cases}\]

dcustom <- function (x) {
  sapply(X = x, 
         FUN = function(x) {
           ifelse(0 < x & x < 1, 
                  (gamma(2+5)/(gamma(2)*gamma(5)))*(x^(2-1)*(1-x)^(5-1)), 
                  0)
         })
}

We first generate a bunch of deviates and then plot their density in teal. Then the true density is overlayed as a yellow line.

df_custom_samples <- 
  data.frame(b = projectq2a(n = 10000, 
                       pdf = dcustom, 
                       a = 0, 
                       b = 1, 
                       C = 2.5))

ggplot(df_custom_samples, aes(x = b)) + 
  geom_density(color = NA, fill = "#35a09c") + 
  stat_function(aes(x = 0), fun = dcustom, color = "#fff735", size = 1) + 
  labs(x = "X", 
       y = "Density") + 
  theme_minimal()

2D Uniform on unit square

We first demo sampling from a uniform distribution defined on the unit square. In order to do so, we must first define a joint probability density function which takes x and y as arguments for the two variables.

\[\text{pdf = }\begin{cases}\frac{1}{(b - a)^2} \text{for } x, y \in [a, b]\\0 \text{ otherwise}\end{cases}\]

d2dunif <- function(x, y, min = 0, max = 1) {
  if(min <= x && x <= max && min <= y && y <= max)
    (max - min)^(-2)
  else
    0
}

Then we could examine what the distribution looks like in 3D by computing the probability densities for each x and y pair in a systematic sampling across the support.

l <- list(x = seq(0, 1, 0.5), 
          y = seq(0, 1, 0.5))

z <- matrix(rep(NA, length(l$x) * length(l$y)), 
            nrow = length(l$x))
for(r in 1:nrow(z)) {
  for(c in 1:ncol(z)) {
    z[r, c] <- d2dunif(x = l$x[r], y = l$y[c])
  }
}
l$z <- z

plot_ly(x = l$x, y = l$y, z = l$z) %>% 
  add_surface() %>% 
  layout(scene = list(xaxis = list(title = "X"), 
                      yaxis = list(title = "Y"), 
                      zaxis = list(title = "Density")))

Finally, we can sample from the joint distribution utilizing rejection sampling and estimate the kernel density.

df <- projectq3a(n = 10000, jpdf = d2dunif, a = 0, b = 1, C = 1)
kd <- with(df, kde2d(x, y, n = 100))
plot_ly(x = kd$x, y = kd$y, z = kd$z) %>% 
  add_surface() %>% 
  layout(scene = list(xaxis = list(title = "X"), 
                      yaxis = list(title = "Y"), 
                      zaxis = list(title = "Density")))

2D Beta on unit square

We first demo sampling from a beta distribution defined on the unit square. In order to do so, we must first define a joint probability density function which takes x and y as arguments for the two variables.

\[\text{pdf = }\begin{cases}\frac{x^{\alpha - 1}(1 - x)^{\beta - 1} + y^{\alpha - 1}(1 - y)^{\beta - 1}}{2B(\alpha, \beta)}\text{ for }0 \le x, y \le 1\\0\text{ otherwise}\end{cases}\]

\[\text{where }B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}, \text{ with }\alpha, \beta \gt 0\]

d2dbeta <- function(x, y, shape1, shape2) {
  if(0 <= x && x <= 1 && 0 <= y && y <= 1)
    (dbeta(x = x, shape1 = shape1, shape2 = shape2) + 
       dbeta(x = y, shape1 = shape1, shape2 = shape2)) / 2
  else
    0
}

Then we could examine what the distribution looks like in 3D by computing the probability densities for each x and y pair in a systematic sampling across the support.

l <- list(x = seq(0, 1, 0.01), 
          y = seq(0, 1, 0.01))

z <- matrix(rep(NA, length(l$x) * length(l$y)), 
            nrow = length(l$x))
for(r in 1:nrow(z)) {
  for(c in 1:ncol(z)) {
    z[r, c] <- d2dbeta(x = l$x[r], y = l$y[c], shape1 = 2, shape2 = 2)
  }
}
l$z <- z

plot_ly(x = l$x, y = l$y, z = l$z) %>% 
  add_surface() %>% 
  layout(scene = list(xaxis = list(title = "X"), 
                      yaxis = list(title = "Y"), 
                      zaxis = list(title = "Density")))

Finally, we can sample from the joint distribution utilizing rejection sampling and estimate the kernel density.

df <- projectq3a(n = 10000, jpdf = d2dbeta, a = 0, b = 1, C = 1.5, shape1 = 2, shape2 = 2)
kd <- with(df, kde2d(x, y, n = 100))
plot_ly(x = kd$x, y = kd$y, z = kd$z) %>% 
  add_surface() %>% 
  layout(scene = list(xaxis = list(title = "X"), 
                      yaxis = list(title = "Y"), 
                      zaxis = list(title = "Density")))

2D Custom on a square

We next demo sampling from a custom joint distribution defined on the square from -1 to +1. Again, we first define the joint pdf function which takes x and y as quantile arguments.

\[\text{pdf = }\begin{cases}\frac{3}{8}(x^2 + y^2) \text{for } x, y \in [-1, +1]\\0 \text{ otherwise}\end{cases}\]

d2dcirclecontour <- function(x, y) {
  if(-1 <= x && x <= 1 && -1 <= y && y <= 1)
    (3/8)*(x^2 + y^2)
  else
    0
}

Then we could examine what the distribution looks like in 3D by computing the probability densities for each x and y pair in a systematic sampling across the support.

l <- list(x = seq(-1, 1, 0.1), 
          y = seq(-1, 1, 0.1))

z <- matrix(rep(NA, length(l$x) * length(l$y)), 
            nrow = length(l$x))
for(r in 1:nrow(z)) {
  for(c in 1:ncol(z)) {
    z[r, c] <- d2dcirclecontour(x = l$x[r], y = l$y[c])
  }
}
l$z <- z

plot_ly(x = l$x, y = l$y, z = l$z) %>% 
  add_surface() %>% 
  layout(scene = list(xaxis = list(title = "X"), 
                      yaxis = list(title = "Y"), 
                      zaxis = list(title = "Density")))

Finally, we can sample from the joint distribution utilizing rejection sampling and estimate the kernel density.

df <- projectq3a(n = 10000, jpdf = d2dcirclecontour, a = -1, b = 1, C = 0.75)
kd <- with(df, kde2d(x, y, n = 100))
plot_ly(x = kd$x, y = kd$y, z = kd$z) %>% 
  add_surface() %>% 
  layout(scene = list(xaxis = list(title = "X"), 
                      yaxis = list(title = "Y"), 
                      zaxis = list(title = "Density")))